program match2sets;
Const
   map = 10; {parameters}
   mp  = 1; {??}
   ndatap  = 3000;
Type
   myreal = real;
   {GAUSSJ, LFIT, MRQMIN, MRQCOF, COVRST...}
   glmma    = ARRAY [1..map] OF REAL;
   glnpbynp   = ARRAY [1..map,1..map] OF REAL;
   {glnpbymp  = ARRAY [1..map,1..mp] OF REAL;}
   glnpbymp   = glnpbynp;
   glnp	      = ARRAY [1..map] OF integer;
   gllista    = glnp;
   glcovar    = glnpbynp;
   glndata    = ARRAY [1..ndatap] OF REAL;
   glncabynca = glnpbymp;
   glnalbynal = glnpbynp;
{MRQMIN}
Var
   covar,alpha			: glcovar;
   ppar,ppard,ppar0,alim0,alim1	: glmma;
   glochisq : REAL;
   glbeta   : glmma;
{MYVAR}
   rf,rw,ala,alao,chi:myreal;
   fo,fo2,fc,fos,fos2,foss,q,qp:glndata;
   dwp,ndat,i,j,k,o,iter,ma,mfit:integer;
   plist:gllista;
   pp:array[1..ndatap] of integer;
   prn:boolean;
   txt:text;
   
procedure funcs_mrq(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
Var
   fc,fcf : myreal;
begin
   fc:=fo2[po]*exp(-a[2]/4*q[po]*q[po]-a[3]/4*qp[po]*qp[po]);
   yfit:=a[1]*fc;
   dyda[1]:=fc;
   dyda[2]:=-yfit*q[po]*q[po]/4;
   dyda[3]:=-yfit*qp[po]*qp[po]/4;
   if (prn) then writeln(po:1,' ',yfit:1:4,' ',fo[po]:1:4,' ',fo2[po]:1:4,' ',a[1]:1:3,' ',a[2]:1:3,' ',a[3]:1:3);
end; { funcs_mrq }

PROCEDURE gaussj(VAR a : glnpbynp; n,np: INTEGER;
		 VAR b : glnpbymp; m,mp: INTEGER);
{ Programs using GAUSSJ must define the types
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnpbymp = ARRAY [1..np,1..mp] OF real;
   glnp = ARRAY [1..np] OF INTEGER;
in the main routine. }
VAR
   big,dum,pivinv: REAL;
   i,icol,irow,j,k,l,ll: INTEGER;
   indxc,indxr,ipiv: glnp;
BEGIN
   FOR j := 1 to n DO BEGIN
      ipiv[j] := 0
   END;
   FOR i := 1 to n DO BEGIN
      big := 0.0;
      FOR j := 1 to n DO BEGIN
         IF (ipiv[j] <> 1) THEN BEGIN
            FOR k := 1 to n DO BEGIN
               IF (ipiv[k] = 0) THEN BEGIN
                  IF (abs(a[j,k]) >= big) THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               END ELSE IF (ipiv[k] > 1) THEN BEGIN
                  {writeln('pause 1 in GAUSSJ - singular matrix'); readln}
               END
            END
         END
      END;
      ipiv[icol] := ipiv[icol]+1;
      IF (irow <> icol) THEN BEGIN
         FOR l := 1 to n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 to m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      END;
      indxr[i] := irow;
      indxc[i] := icol;
      IF (a[icol,icol] = 0.0) THEN BEGIN
	 writeln('pause 2 in GAUSSJ - singular matrix'); {readln}
	 {writeln('pause 2 in GAUSSJ - singular matrix : ma=',plist[icol]:1,' pdw psca ',pdw:1,' ',psca:1);
	  if (plist[icol]<(psca-1)) then 
	  writeln('ORB ',p2o[plist[icol],1]:1,' CHEM ',p2o[plist[icol],2]:1) else
	  writeln('zero=',p2o[plist[icol],1]:1,'[1=sca 2=pdw 3=aconst] :  ',p2o[plist[icol],2]:1);
	  }
	 halt;
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 to n DO BEGIN
         a[icol,l] := a[icol,l]*pivinv
      END;
      FOR l := 1 to m DO BEGIN
         b[icol,l] := b[icol,l]*pivinv
      END;
      FOR ll := 1 to n DO BEGIN
         IF (ll <> icol) THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 to n DO BEGIN
               a[ll,l] := a[ll,l]-a[icol,l]*dum
            END;
            FOR l := 1 to m DO BEGIN
               b[ll,l] := b[ll,l]-b[icol,l]*dum
            END
         END
      END
   END;
   FOR l := n DOWNTO 1 DO BEGIN
      IF (indxr[l] <> indxc[l]) THEN BEGIN
         FOR k := 1 to n DO BEGIN
            dum := a[k,indxr[l]];
            a[k,indxr[l]] := a[k,indxc[l]];
            a[k,indxc[l]] := dum
         END
      END
   END
END;
PROCEDURE covsrt(VAR covar: glcovar; ncvm: INTEGER; ma: INTEGER;
       lista: gllista; mfit: INTEGER);
{ Programs using routine COVSRT must define the types
TYPE
   glcovar = ARRAY [1..ncvm,1..ncvm] OF real;
   gllista = ARRAY [1..mfit] OF INTEGER;
in the calling program. }
VAR
   j,i: INTEGER;
   swap: REAL;
BEGIN
   FOR j := 1 to ma-1 DO BEGIN
      FOR i := j+1 to ma DO BEGIN
         covar[i,j] := 0.0
      END
   END;
   FOR i := 1 to mfit-1 DO BEGIN
      FOR j := i+1 to mfit DO BEGIN
         IF (lista[j] > lista[i])  THEN BEGIN
            covar[lista[j],lista[i]] := covar[i,j]
         END ELSE BEGIN
            covar[lista[i],lista[j]] := covar[i,j]
         END
      END
   END;
   swap := covar[1,1];
   FOR j := 1 to ma DO BEGIN
      covar[1,j] := covar[j,j];
      covar[j,j] := 0.0
   END;
   covar[lista[1],lista[1]] := swap;
   FOR j := 2 to mfit DO BEGIN
      covar[lista[j],lista[j]] := covar[1,j]
   END;
   FOR j := 2 to ma DO BEGIN
      FOR i := 1 to j-1 DO BEGIN
         covar[i,j] := covar[j,i]
      END
   END
END;
PROCEDURE mrqcof(x,y,sig: glndata; ndata: INTEGER;
       VAR a: glmma; mma: INTEGER; lista: gllista;
       mfit: INTEGER; VAR alpha: glnalbynal;
       VAR beta: glmma; nalp: INTEGER; VAR chisq: REAL);
VAR
   k,j,i	    : INTEGER;
   ymod,wt,sig2i,dy : REAL;
   dyda		    : glmma;
   da0,da1	    : real;
BEGIN
   FOR j := 1 to mfit DO BEGIN
      FOR k := 1 to j DO BEGIN
         alpha[j,k] := 0.0
      END;
      beta[j] := 0.0
   END;
   chisq := 0.0;
   {writeln('mrqcof ',mma:1);}
   FOR i := 1 to ndata DO BEGIN
      funcs_mrq(i,a,ymod,dyda,mma);
      sig2i := 1.0/(sig[i]*sig[i]);
      dy := y[i]-ymod;
      FOR j := 1 to mfit DO BEGIN
         wt := dyda[lista[j]]*sig2i;
         FOR k := 1 to j DO BEGIN
            alpha[j,k] := alpha[j,k]+wt*dyda[lista[k]]
         END;
         beta[j] := beta[j]+dy*wt
      END;
      chisq := chisq+dy*dy*sig2i
   END;
   {writeln('CHISQ_TRIAL: ',chisq/ndata:1:5);}
   FOR j := 2 to mfit DO BEGIN
      FOR k := 1 to j-1 DO BEGIN
         alpha[k,j] := alpha[j,k]
      END
   END
END;
PROCEDURE mrqmin(x,y,sig : glndata; ndata: INTEGER;VAR a: glmma; mma: INTEGER; lista: gllista; mfit: INTEGER; VAR covar,alpha: glncabynca;nca: INTEGER; VAR chisq,alamda: REAL);
VAR
   k,kk,j,ihit : INTEGER;
   atry,da     : glmma; oneda: glncabynca;
   goto99      : boolean;
BEGIN
   IF (alamda < 0.0) THEN BEGIN
      kk := mfit+1;
      FOR j := 1 TO mma DO BEGIN
	 ihit := 0;
	 FOR k := 1 TO mfit DO BEGIN
	    IF (lista[k] = j) THEN ihit := ihit+1
	 END;
	 IF (ihit = 0) THEN BEGIN
	    lista[kk] := j;  kk := kk+1
	 END ELSE IF (ihit > 1) THEN BEGIN
	       writeln('pause 1 in routine MRQMIN');
	     writeln('Improper permutation in LISTA'); readln
	 END
      END;
      IF (kk <> (mma+1)) THEN BEGIN
	 writeln('pause 2 in routine MRQMIN');
	 writeln('Improper permutation in LISTA'); readln
      END;
      alamda := 0.001;
      mrqcof(x,y,sig,ndata,a,mma,lista,mfit,alpha,glbeta,nca,chisq);
      glochisq := chisq;
      FOR j := 1 TO mma DO BEGIN
	   atry[j] := a[j]
      END
   END;
   FOR j := 1 TO mfit DO BEGIN
      FOR k := 1 TO mfit DO BEGIN
	 covar[j,k] := alpha[j,k]
      END;
      covar[j,j] := alpha[j,j]*(1.0+alamda);
      oneda[j,1] := glbeta[j]
   END;
   gaussj(covar,mfit,nca,oneda,1,1);
   FOR j := 1 TO mfit DO da[j] := oneda[j,1];
   goto99:=false;
   IF (alamda < 0.0) THEN BEGIN
      covsrt(covar,nca,mma,lista,mfit);
      GOTO99:=true;
   END;

   if not(goto99) then
   begin      
      FOR j := 1 TO mfit DO BEGIN
	 atry[lista[j]] := a[lista[j]]+da[j]
      END;

      mrqcof(x,y,sig,ndata,atry,mma,lista,mfit,covar,da,nca,chisq);
      IF (chisq < glochisq) THEN BEGIN
	 alamda := 0.1 * alamda;  glochisq := chisq;
	 FOR j := 1 TO mfit DO BEGIN
	    FOR k := 1 TO mfit DO BEGIN
	       alpha[j,k] := covar[j,k]
	    END;
	    glbeta[j] := da[j];
	    a[lista[j]] := atry[lista[j]]
	 END
      END
      ELSE BEGIN
	 {for j:=1 to mfit do stn[lista[j]]:=-1;}
	 alamda := 10.0*alamda;  chisq := glochisq
      END;
   end;
END; { mrqmin }
procedure rwfac(fo,fc,sig : glndata;ndata:INTEGER; Var rfac,rwfac,mchi:myreal);
var
   i,j			   : INTEGER;
   dy,sum,dif,sumw,difw,wt : MYREAL;
BEGIN
{ 
RELIABILITY FACTORS (f=|F|):
R = Sum(|f_o-f_c|)/Sum(f_o)
wR = Sum(w*|f_o-f_c|)/Sum(w*f_o)
w =1/sig^2
}
   dif:=0; sum:=0;
   difw:=0; sumw:=0;
   mchi:=0;
   for i:=1 to ndata do
   begin
      sum:=sum+fo[i];
      dy:=abs(fo[i]-fc[i]);
      dif:=dif+dy;
      wt:=1/sig[i]/sig[i];
      sumw:=sumw+wt*fo[i];
      difw:=difw+wt*abs(fo[i]-fc[i]);
      mchi:=mchi+dy*dy*wt;
   end;
   rfac:=dif/sum;
   rwfac:=difw/sumw;
   {mchi:=mchi/ndata;}
   mchi:=mchi/ndat;
END; { rwfac }
BEGIN
   i:=0;
   repeat
      i:=i+1;
      readln(fo[i],fo2[i],fos[i],fos2[i],q[i],qp[i],o,pp[i]);
      foss[i]:=fos2[i];
      {fo[i]:=sqrt(fo[i]);}
      {writeln(i:1,' ',fo[i]:1:4,' ',fos[i]:1:4);}
   until eof;
   ndat:=i-1;
   {writeln(ndat:1);}
   {set parameters}
   ppar[1]:=1; {scale}
   ppar[2]:=0; {dw}
   ppar[3]:=0; {perpdw}
   ppar[4]:=0;
   ma:=3;
   for o:=1 to ma do plist[o]:=o;
   mfit:=3;
   iter:=0; ala:=-0.1; alao:=ala;
   prn:=false;
   repeat
      iter:=iter+1;
      mrqmin(q,fo,foss,ndat,ppar,ma,plist,ma,covar,alpha,ma,chi,ala);
   until (iter<=200);
   write('chi=',chi/ndat:1:3,'  ppar: ');
   for o:=1 to mfit do write(ppar[o]:1:5,' '); writeln;
   prn:=false;
   for i:=1 to ndat do begin
      funcs_mrq(i,ppar,fc[i],ppard,4);
   end;
   rwfac(fo,fc,fos,ndat,rf,rw,chi);
   writeln('wR R chi : ',rw:1:3,' ',rf:1:3,' ',chi:1:3);
   rewrite(txt,'cdata.rep');
   for i:=1 to ndat do begin
      writeln(txt,fo[i]:1:6,' ',fc[i]:1:6,' ',fos[i]:1:6,' ',fos2[i]:1:6,' ',q[i]:1:4,' ',qp[i]:1:4,' ',pp[i]:1);
   end;
   close(txt);
END.